home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1993 July / InfoMagic USENET CD-ROM July 1993.ISO / sources / misc / volume24 / gnucalc / part21 < prev    next >
Encoding:
Text File  |  1991-10-31  |  55.2 KB  |  1,868 lines

  1. Newsgroups: comp.sources.misc
  2. From: daveg@synaptics.com (David Gillespie)
  3. Subject:  v24i069:  gnucalc - GNU Emacs Calculator, v2.00, Part21/56
  4. Message-ID: <1991Oct31.072644.17973@sparky.imd.sterling.com>
  5. X-Md4-Signature: 935912531147058885be56ee6ad9de74
  6. Date: Thu, 31 Oct 1991 07:26:44 GMT
  7. Approved: kent@sparky.imd.sterling.com
  8.  
  9. Submitted-by: daveg@synaptics.com (David Gillespie)
  10. Posting-number: Volume 24, Issue 69
  11. Archive-name: gnucalc/part21
  12. Environment: Emacs
  13. Supersedes: gmcalc: Volume 13, Issue 27-45
  14.  
  15. ---- Cut Here and unpack ----
  16. #!/bin/sh
  17. # do not concatenate these parts, unpack them in order with /bin/sh
  18. # file calc-mat.el continued
  19. #
  20. if test ! -r _shar_seq_.tmp; then
  21.     echo 'Please unpack part 1 first!'
  22.     exit 1
  23. fi
  24. (read Scheck
  25.  if test "$Scheck" != 21; then
  26.     echo Please unpack part "$Scheck" next!
  27.     exit 1
  28.  else
  29.     exit 0
  30.  fi
  31. ) < _shar_seq_.tmp || exit 1
  32. if test ! -f _shar_wnt_.tmp; then
  33.     echo 'x - still skipping calc-mat.el'
  34. else
  35. echo 'x - continuing file calc-mat.el'
  36. sed 's/^X//' << 'SHAR_EOF' >> 'calc-mat.el' &&
  37. X                        (nth j (nth k lu))))
  38. X        k (1+ k)))
  39. X    (setcar (nthcdr j (nth i lu)) sum)
  40. X    (let ((dum (math-abs-approx sum)))
  41. X      (if (Math-lessp big dum)
  42. X          (setq big dum
  43. X            imax i)))
  44. X    (setq i (1+ i)))
  45. X      (if (> imax j)
  46. X      (setq lu (math-swap-rows lu j imax)
  47. X        d (- d)))
  48. X      (setq index (cons imax index))
  49. X      (let ((pivot (nth j (nth j lu))))
  50. X    (if (math-zerop pivot)
  51. X        (throw 'singular nil)
  52. X      (setq i j)
  53. X      (while (<= (setq i (1+ i)) n)
  54. X        (setcar (nthcdr j (nth i lu))
  55. X            (math-div (nth j (nth i lu)) pivot)))))
  56. X      (setq j (1+ j)))
  57. X    (list lu (nreverse index) d))
  58. )
  59. X
  60. (defun math-swap-rows (m r1 r2)
  61. X  (or (= r1 r2)
  62. X      (let* ((r1prev (nthcdr (1- r1) m))
  63. X         (row1 (cdr r1prev))
  64. X         (r2prev (nthcdr (1- r2) m))
  65. X         (row2 (cdr r2prev))
  66. X         (r2next (cdr row2)))
  67. X    (setcdr r2prev row1)
  68. X    (setcdr r1prev row2)
  69. X    (setcdr row2 (cdr row1))
  70. X    (setcdr row1 r2next)))
  71. X  m
  72. )
  73. X
  74. X
  75. (defun math-lud-solve (lud b &optional need)
  76. X  (if lud
  77. X      (let* ((x (math-copy-matrix b))
  78. X         (n (1- (length x)))
  79. X         (m (1- (length (nth 1 x))))
  80. X         (lu (car lud))
  81. X         (col 1)
  82. X         i j ip ii index sum)
  83. X    (while (<= col m)
  84. X      (math-working "LUD solver step" col)
  85. X      (setq i 1
  86. X        ii nil
  87. X        index (nth 1 lud))
  88. X      (while (<= i n)
  89. X        (setq ip (car index)
  90. X          index (cdr index)
  91. X          sum (nth col (nth ip x)))
  92. X        (setcar (nthcdr col (nth ip x)) (nth col (nth i x)))
  93. X        (if (null ii)
  94. X        (or (math-zerop sum)
  95. X            (setq ii i))
  96. X          (setq j ii)
  97. X          (while (< j i)
  98. X        (setq sum (math-sub sum (math-mul (nth j (nth i lu))
  99. X                          (nth col (nth j x))))
  100. X              j (1+ j))))
  101. X        (setcar (nthcdr col (nth i x)) sum)
  102. X        (setq i (1+ i)))
  103. X      (while (>= (setq i (1- i)) 1)
  104. X        (setq sum (nth col (nth i x))
  105. X          j i)
  106. X        (while (<= (setq j (1+ j)) n)
  107. X          (setq sum (math-sub sum (math-mul (nth j (nth i lu))
  108. X                        (nth col (nth j x))))))
  109. X        (setcar (nthcdr col (nth i x))
  110. X            (math-div sum (nth i (nth i lu)))))
  111. X      (setq col (1+ col)))
  112. X    x)
  113. X    (and need
  114. X     (math-reject-arg need "*Singular matrix")))
  115. )
  116. X
  117. (defun calcFunc-lud (m)
  118. X  (if (math-square-matrixp m)
  119. X      (or (math-with-extra-prec 2
  120. X        (let ((lud (math-matrix-lud m)))
  121. X          (and lud
  122. X           (let* ((lmat (math-copy-matrix (car lud)))
  123. X              (umat (math-copy-matrix (car lud)))
  124. X              (n (1- (length (car lud))))
  125. X              (perm (calcFunc-idn 1 n))
  126. X              i (j 1))
  127. X             (while (<= j n)
  128. X               (setq i 1)
  129. X               (while (< i j)
  130. X             (setcar (nthcdr j (nth i lmat)) 0)
  131. X             (setq i (1+ i)))
  132. X               (setcar (nthcdr j (nth j lmat)) 1)
  133. X               (while (<= (setq i (1+ i)) n)
  134. X             (setcar (nthcdr j (nth i umat)) 0))
  135. X               (setq j (1+ j)))
  136. X             (while (>= (setq j (1- j)) 1)
  137. X               (let ((pos (nth (1- j) (nth 1 lud))))
  138. X             (or (= pos j)
  139. X                 (setq perm (math-swap-rows perm j pos)))))
  140. X             (list 'vec perm lmat umat)))))
  141. X      (math-reject-arg m "*Singular matrix"))
  142. X    (math-reject-arg m 'square-matrixp))
  143. )
  144. X
  145. SHAR_EOF
  146. echo 'File calc-mat.el is complete' &&
  147. chmod 0644 calc-mat.el ||
  148. echo 'restore of calc-mat.el failed'
  149. Wc_c="`wc -c < 'calc-mat.el'`"
  150. test 10372 -eq "$Wc_c" ||
  151.     echo 'calc-mat.el: original size 10372, current size' "$Wc_c"
  152. rm -f _shar_wnt_.tmp
  153. fi
  154. # ============= calc-math.el ==============
  155. if test -f 'calc-math.el' -a X"$1" != X"-c"; then
  156.     echo 'x - skipping calc-math.el (File already exists)'
  157.     rm -f _shar_wnt_.tmp
  158. else
  159. > _shar_wnt_.tmp
  160. echo 'x - extracting calc-math.el (Text)'
  161. sed 's/^X//' << 'SHAR_EOF' > 'calc-math.el' &&
  162. ;; Calculator for GNU Emacs, part II [calc-math.el]
  163. ;; Copyright (C) 1990, 1991 Free Software Foundation, Inc.
  164. ;; Written by Dave Gillespie, daveg@csvax.cs.caltech.edu.
  165. X
  166. ;; This file is part of GNU Emacs.
  167. X
  168. ;; GNU Emacs is distributed in the hope that it will be useful,
  169. ;; but WITHOUT ANY WARRANTY.  No author or distributor
  170. ;; accepts responsibility to anyone for the consequences of using it
  171. ;; or for whether it serves any particular purpose or works at all,
  172. ;; unless he says so in writing.  Refer to the GNU Emacs General Public
  173. ;; License for full details.
  174. X
  175. ;; Everyone is granted permission to copy, modify and redistribute
  176. ;; GNU Emacs, but only under the conditions described in the
  177. ;; GNU Emacs General Public License.   A copy of this license is
  178. ;; supposed to have been given to you along with GNU Emacs so you
  179. ;; can know your rights and responsibilities.  It should be in a
  180. ;; file named COPYING.  Among other things, the copyright notice
  181. ;; and this notice must be preserved on all copies.
  182. X
  183. X
  184. X
  185. ;; This file is autoloaded from calc-ext.el.
  186. (require 'calc-ext)
  187. X
  188. (require 'calc-macs)
  189. X
  190. (defun calc-Need-calc-math () nil)
  191. X
  192. X
  193. (defun calc-sqrt (arg)
  194. X  (interactive "P")
  195. X  (calc-slow-wrapper
  196. X   (if (calc-is-inverse)
  197. X       (calc-unary-op "^2" 'calcFunc-sqr arg)
  198. X     (calc-unary-op "sqrt" 'calcFunc-sqrt arg)))
  199. )
  200. X
  201. (defun calc-isqrt (arg)
  202. X  (interactive "P")
  203. X  (calc-slow-wrapper
  204. X   (if (calc-is-inverse)
  205. X       (calc-unary-op "^2" 'calcFunc-sqr arg)
  206. X     (calc-unary-op "isqt" 'calcFunc-isqrt arg)))
  207. )
  208. X
  209. X
  210. (defun calc-hypot (arg)
  211. X  (interactive "P")
  212. X  (calc-slow-wrapper
  213. X   (calc-binary-op "hypt" 'calcFunc-hypot arg))
  214. )
  215. X
  216. (defun calc-ln (arg)
  217. X  (interactive "P")
  218. X  (calc-invert-func)
  219. X  (calc-exp arg)
  220. )
  221. X
  222. (defun calc-log10 (arg)
  223. X  (interactive "P")
  224. X  (calc-hyperbolic-func)
  225. X  (calc-ln arg)
  226. )
  227. X
  228. (defun calc-log (arg)
  229. X  (interactive "P")
  230. X  (calc-slow-wrapper
  231. X   (if (calc-is-inverse)
  232. X       (calc-binary-op "alog" 'calcFunc-alog arg)
  233. X     (calc-binary-op "log" 'calcFunc-log arg)))
  234. )
  235. X
  236. (defun calc-ilog (arg)
  237. X  (interactive "P")
  238. X  (calc-slow-wrapper
  239. X   (if (calc-is-inverse)
  240. X       (calc-binary-op "alog" 'calcFunc-alog arg)
  241. X     (calc-binary-op "ilog" 'calcFunc-ilog arg)))
  242. )
  243. X
  244. (defun calc-lnp1 (arg)
  245. X  (interactive "P")
  246. X  (calc-invert-func)
  247. X  (calc-expm1 arg)
  248. )
  249. X
  250. (defun calc-exp (arg)
  251. X  (interactive "P")
  252. X  (calc-slow-wrapper
  253. X   (if (calc-is-hyperbolic)
  254. X       (if (calc-is-inverse)
  255. X       (calc-unary-op "lg10" 'calcFunc-log10 arg)
  256. X     (calc-unary-op "10^" 'calcFunc-exp10 arg))
  257. X     (if (calc-is-inverse)
  258. X     (calc-unary-op "ln" 'calcFunc-ln arg)
  259. X       (calc-unary-op "exp" 'calcFunc-exp arg))))
  260. )
  261. X
  262. (defun calc-expm1 (arg)
  263. X  (interactive "P")
  264. X  (calc-slow-wrapper
  265. X   (if (calc-is-inverse)
  266. X       (calc-unary-op "ln+1" 'calcFunc-lnp1 arg)
  267. X     (calc-unary-op "ex-1" 'calcFunc-expm1 arg)))
  268. )
  269. X
  270. (defun calc-pi ()
  271. X  (interactive)
  272. X  (calc-slow-wrapper
  273. X   (if (calc-is-inverse)
  274. X       (if (calc-is-hyperbolic)
  275. X       (if calc-symbolic-mode
  276. X           (calc-pop-push-record 0 "phi" '(var phi var-phi))
  277. X         (calc-pop-push-record 0 "phi" (math-phi)))
  278. X     (if calc-symbolic-mode
  279. X         (calc-pop-push-record 0 "gmma" '(var gamma var-gamma))
  280. X       (calc-pop-push-record 0 "gmma" (math-gamma-const))))
  281. X     (if (calc-is-hyperbolic)
  282. X     (if calc-symbolic-mode
  283. X         (calc-pop-push-record 0 "e" '(var e var-e))
  284. X       (calc-pop-push-record 0 "e" (math-e)))
  285. X       (if calc-symbolic-mode
  286. X       (calc-pop-push-record 0 "pi" '(var pi var-pi))
  287. X     (calc-pop-push-record 0 "pi" (math-pi))))))
  288. )
  289. X
  290. (defun calc-sin (arg)
  291. X  (interactive "P")
  292. X  (calc-slow-wrapper
  293. X   (if (calc-is-hyperbolic)
  294. X       (if (calc-is-inverse)
  295. X       (calc-unary-op "asnh" 'calcFunc-arcsinh arg)
  296. X     (calc-unary-op "sinh" 'calcFunc-sinh arg))
  297. X     (if (calc-is-inverse)
  298. X     (calc-unary-op "asin" 'calcFunc-arcsin arg)
  299. X       (calc-unary-op "sin" 'calcFunc-sin arg))))
  300. )
  301. X
  302. (defun calc-arcsin (arg)
  303. X  (interactive "P")
  304. X  (calc-invert-func)
  305. X  (calc-sin arg)
  306. )
  307. X
  308. (defun calc-sinh (arg)
  309. X  (interactive "P")
  310. X  (calc-hyperbolic-func)
  311. X  (calc-sin arg)
  312. )
  313. X
  314. (defun calc-arcsinh (arg)
  315. X  (interactive "P")
  316. X  (calc-invert-func)
  317. X  (calc-hyperbolic-func)
  318. X  (calc-sin arg)
  319. )
  320. X
  321. (defun calc-cos (arg)
  322. X  (interactive "P")
  323. X  (calc-slow-wrapper
  324. X   (if (calc-is-hyperbolic)
  325. X       (if (calc-is-inverse)
  326. X       (calc-unary-op "acsh" 'calcFunc-arccosh arg)
  327. X     (calc-unary-op "cosh" 'calcFunc-cosh arg))
  328. X     (if (calc-is-inverse)
  329. X     (calc-unary-op "acos" 'calcFunc-arccos arg)
  330. X       (calc-unary-op "cos" 'calcFunc-cos arg))))
  331. )
  332. X
  333. (defun calc-arccos (arg)
  334. X  (interactive "P")
  335. X  (calc-invert-func)
  336. X  (calc-cos arg)
  337. )
  338. X
  339. (defun calc-cosh (arg)
  340. X  (interactive "P")
  341. X  (calc-hyperbolic-func)
  342. X  (calc-cos arg)
  343. )
  344. X
  345. (defun calc-arccosh (arg)
  346. X  (interactive "P")
  347. X  (calc-invert-func)
  348. X  (calc-hyperbolic-func)
  349. X  (calc-cos arg)
  350. )
  351. X
  352. (defun calc-sincos ()
  353. X  (interactive)
  354. X  (calc-slow-wrapper
  355. X   (if (calc-is-inverse)
  356. X       (calc-enter-result 1 "asnc" (list 'calcFunc-arcsincos (calc-top-n 1)))
  357. X     (calc-enter-result 1 "sncs" (list 'calcFunc-sincos (calc-top-n 1)))))
  358. )
  359. X
  360. (defun calc-tan (arg)
  361. X  (interactive "P")
  362. X  (calc-slow-wrapper
  363. X   (if (calc-is-hyperbolic)
  364. X       (if (calc-is-inverse)
  365. X       (calc-unary-op "atnh" 'calcFunc-arctanh arg)
  366. X     (calc-unary-op "tanh" 'calcFunc-tanh arg))
  367. X     (if (calc-is-inverse)
  368. X     (calc-unary-op "atan" 'calcFunc-arctan arg)
  369. X       (calc-unary-op "tan" 'calcFunc-tan arg))))
  370. )
  371. X
  372. (defun calc-arctan (arg)
  373. X  (interactive "P")
  374. X  (calc-invert-func)
  375. X  (calc-tan arg)
  376. )
  377. X
  378. (defun calc-tanh (arg)
  379. X  (interactive "P")
  380. X  (calc-hyperbolic-func)
  381. X  (calc-tan arg)
  382. )
  383. X
  384. (defun calc-arctanh (arg)
  385. X  (interactive "P")
  386. X  (calc-invert-func)
  387. X  (calc-hyperbolic-func)
  388. X  (calc-tan arg)
  389. )
  390. X
  391. (defun calc-arctan2 ()
  392. X  (interactive)
  393. X  (calc-slow-wrapper
  394. X   (calc-enter-result 2 "atn2" (cons 'calcFunc-arctan2 (calc-top-list-n 2))))
  395. )
  396. X
  397. (defun calc-conj (arg)
  398. X  (interactive "P")
  399. X  (calc-wrapper
  400. X   (calc-unary-op "conj" 'calcFunc-conj arg))
  401. )
  402. X
  403. (defun calc-imaginary ()
  404. X  (interactive)
  405. X  (calc-slow-wrapper
  406. X   (calc-pop-push-record 1 "i*" (math-imaginary (calc-top-n 1))))
  407. )
  408. X
  409. X
  410. X
  411. (defun calc-to-degrees (arg)
  412. X  (interactive "P")
  413. X  (calc-wrapper
  414. X   (calc-unary-op ">deg" 'calcFunc-deg arg))
  415. )
  416. X
  417. (defun calc-to-radians (arg)
  418. X  (interactive "P")
  419. X  (calc-wrapper
  420. X   (calc-unary-op ">rad" 'calcFunc-rad arg))
  421. )
  422. X
  423. X
  424. (defun calc-degrees-mode ()
  425. X  (interactive)
  426. X  (calc-wrapper
  427. X   (calc-change-mode 'calc-angle-mode 'deg)
  428. X   (message "Angles measured in degrees."))
  429. )
  430. X
  431. (defun calc-radians-mode ()
  432. X  (interactive)
  433. X  (calc-wrapper
  434. X   (calc-change-mode 'calc-angle-mode 'rad)
  435. X   (message "Angles measured in radians."))
  436. )
  437. X
  438. X
  439. ;;; Compute the integer square-root floor(sqrt(A)).  A > 0.  [I I] [Public]
  440. ;;; This method takes advantage of the fact that Newton's method starting
  441. ;;; with an overestimate always works, even using truncating integer division!
  442. (defun math-isqrt (a)
  443. X  (cond ((Math-zerop a) a)
  444. X    ((not (math-natnump a))
  445. X     (math-reject-arg a 'natnump))
  446. X    ((integerp a)
  447. X     (math-isqrt-small a))
  448. X    (t
  449. X     (math-normalize (cons 'bigpos (cdr (math-isqrt-bignum (cdr a)))))))
  450. )
  451. X
  452. (defun calcFunc-isqrt (a)
  453. X  (if (math-realp a)
  454. X      (math-isqrt (math-floor a))
  455. X    (math-floor (math-sqrt a)))
  456. )
  457. X
  458. X
  459. ;;; This returns (flag . result) where the flag is T if A is a perfect square.
  460. (defun math-isqrt-bignum (a)   ; [P.l L]
  461. X  (let ((len (length a)))
  462. X    (if (= (% len 2) 0)
  463. X    (let* ((top (nthcdr (- len 2) a)))
  464. X      (math-isqrt-bignum-iter
  465. X       a
  466. X       (math-scale-bignum-3
  467. X        (math-bignum-big
  468. X         (1+ (math-isqrt-small
  469. X          (+ (* (nth 1 top) 1000) (car top)))))
  470. X        (1- (/ len 2)))))
  471. X      (let* ((top (nth (1- len) a)))
  472. X    (math-isqrt-bignum-iter
  473. X     a
  474. X     (math-scale-bignum-3
  475. X      (list (1+ (math-isqrt-small top)))
  476. X      (/ len 2))))))
  477. )
  478. X
  479. (defun math-isqrt-bignum-iter (a guess)   ; [l L l]
  480. X  (math-working "isqrt" (cons 'bigpos guess))
  481. X  (let* ((q (math-div-bignum a guess))
  482. X     (s (math-add-bignum (car q) guess))
  483. X     (g2 (math-div2-bignum s))
  484. X     (comp (math-compare-bignum g2 guess)))
  485. X    (if (< comp 0)
  486. X    (math-isqrt-bignum-iter a g2)
  487. X      (cons (and (= comp 0)
  488. X         (math-zerop-bignum (cdr q))
  489. X         (= (% (car s) 2) 0))
  490. X        guess)))
  491. )
  492. X
  493. (defun math-zerop-bignum (a)
  494. X  (and (eq (car a) 0)
  495. X       (progn
  496. X     (while (eq (car (setq a (cdr a))) 0))
  497. X     (null a)))
  498. )
  499. X
  500. (defun math-scale-bignum-3 (a n)   ; [L L S]
  501. X  (while (> n 0)
  502. X    (setq a (cons 0 a)
  503. X      n (1- n)))
  504. X  a
  505. )
  506. X
  507. (defun math-isqrt-small (a)   ; A > 0.  [S S]
  508. X  (let ((g (cond ((>= a 10000) 1000)
  509. X         ((>= a 100) 100)
  510. X         (t 10)))
  511. X    g2)
  512. X    (while (< (setq g2 (/ (+ g (/ a g)) 2)) g)
  513. X      (setq g g2))
  514. X    g)
  515. )
  516. X
  517. X
  518. X
  519. X
  520. ;;; Compute the square root of a number.
  521. ;;; [T N] if possible, else [F N] if possible, else [C N].  [Public]
  522. (defun math-sqrt (a)
  523. X  (or
  524. X   (and (Math-zerop a) a)
  525. X   (and (math-known-nonposp a)
  526. X    (math-imaginary (math-sqrt (math-neg a))))
  527. X   (and (integerp a)
  528. X    (let ((sqrt (math-isqrt-small a)))
  529. X      (if (= (* sqrt sqrt) a)
  530. X          sqrt
  531. X        (if calc-symbolic-mode
  532. X        (list 'calcFunc-sqrt a)
  533. X          (math-sqrt-float (math-float a) (math-float sqrt))))))
  534. X   (and (eq (car-safe a) 'bigpos)
  535. X    (let* ((res (math-isqrt-bignum (cdr a)))
  536. X           (sqrt (math-normalize (cons 'bigpos (cdr res)))))
  537. X      (if (car res)
  538. X          sqrt
  539. X        (if calc-symbolic-mode
  540. X        (list 'calcFunc-sqrt a)
  541. X          (math-sqrt-float (math-float a) (math-float sqrt))))))
  542. X   (and (eq (car-safe a) 'frac)
  543. X    (let* ((num-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 1 a)))))
  544. X           (num-sqrt (math-normalize (cons 'bigpos (cdr num-res))))
  545. X           (den-res (math-isqrt-bignum (cdr (Math-bignum-test (nth 2 a)))))
  546. X           (den-sqrt (math-normalize (cons 'bigpos (cdr den-res)))))
  547. X      (if (and (car num-res) (car den-res))
  548. X          (list 'frac num-sqrt den-sqrt)
  549. X        (if calc-symbolic-mode
  550. X        (if (or (car num-res) (car den-res))
  551. X            (math-div (if (car num-res)
  552. X                  num-sqrt (list 'calcFunc-sqrt (nth 1 a)))
  553. X                  (if (car den-res)
  554. X                  den-sqrt (list 'calcFunc-sqrt (nth 2 a))))
  555. X          (list 'calcFunc-sqrt a))
  556. X          (math-sqrt-float (math-float a)
  557. X                   (math-div (math-float num-sqrt) den-sqrt))))))
  558. X   (and (eq (car-safe a) 'float)
  559. X    (if calc-symbolic-mode
  560. X        (if (= (% (nth 2 a) 2) 0)
  561. X        (let ((res (math-isqrt-bignum
  562. X                (cdr (Math-bignum-test (nth 1 a))))))
  563. X          (if (car res)
  564. X              (math-make-float (math-normalize
  565. X                    (cons 'bigpos (cdr res)))
  566. X                       (/ (nth 2 a) 2))
  567. X            (signal 'inexact-result nil)))
  568. X          (signal 'inexact-result nil))
  569. X      (math-sqrt-float a)))
  570. X   (and (eq (car-safe a) 'cplx)
  571. X    (math-with-extra-prec 2
  572. X      (let* ((d (math-abs a))
  573. X         (imag (math-sqrt (math-mul (math-sub d (nth 1 a))
  574. X                        '(float 5 -1)))))
  575. X        (list 'cplx
  576. X          (math-sqrt (math-mul (math-add d (nth 1 a)) '(float 5 -1)))
  577. X          (if (math-negp (nth 2 a)) (math-neg imag) imag)))))
  578. X   (and (eq (car-safe a) 'polar)
  579. X    (list 'polar
  580. X          (math-sqrt (nth 1 a))
  581. X          (math-mul (nth 2 a) '(float 5 -1))))
  582. X   (and (eq (car-safe a) 'sdev)
  583. X    (let ((sqrt (math-sqrt (nth 1 a))))
  584. X      (math-make-sdev sqrt
  585. X              (math-div (nth 2 a) (math-mul sqrt 2)))))
  586. X   (and (eq (car-safe a) 'intv)
  587. X    (not (math-negp (nth 2 a)))
  588. X    (math-make-intv (nth 1 a) (math-sqrt (nth 2 a)) (math-sqrt (nth 3 a))))
  589. X   (and (eq (car-safe a) '*)
  590. X    (or (math-known-nonnegp (nth 1 a))
  591. X        (math-known-nonnegp (nth 2 a)))
  592. X    (math-mul (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
  593. X   (and (eq (car-safe a) '/)
  594. X    (or (and (math-known-nonnegp (nth 2 a))
  595. X         (math-div (math-sqrt (nth 1 a)) (math-sqrt (nth 2 a))))
  596. X        (and (math-known-nonnegp (nth 1 a))
  597. X         (not (math-equal-int (nth 1 a) 1))
  598. X         (math-mul (math-sqrt (nth 1 a))
  599. X               (math-sqrt (math-div 1 (nth 2 a)))))))
  600. X   (and (eq (car-safe a) '^)
  601. X    (math-known-evenp (nth 2 a))
  602. X    (math-known-realp (nth 1 a))
  603. X    (math-abs (math-pow (nth 1 a) (math-div (nth 2 a) 2))))
  604. X   (let ((inf (math-infinitep a)))
  605. X     (and inf
  606. X      (math-mul (math-sqrt (math-infinite-dir a inf)) inf)))
  607. X   (progn
  608. X     (calc-record-why 'numberp a)
  609. X     (list 'calcFunc-sqrt a)))
  610. )
  611. (fset 'calcFunc-sqrt (symbol-function 'math-sqrt))
  612. X
  613. (defun math-infinite-dir (a &optional inf)
  614. X  (or inf (setq inf (math-infinitep a)))
  615. X  (math-normalize (math-expr-subst a inf 1))
  616. )
  617. X
  618. (defun math-sqrt-float (a &optional guess)   ; [F F F]
  619. X  (if calc-symbolic-mode
  620. X      (signal 'inexact-result nil)
  621. X    (math-with-extra-prec 1 (math-sqrt-raw a guess)))
  622. )
  623. X
  624. (defun math-sqrt-raw (a &optional guess)   ; [F F F]
  625. X  (if (not (Math-posp a))
  626. X      (math-sqrt a)
  627. X    (if (null guess)
  628. X    (let ((ldiff (- (math-numdigs (nth 1 a)) 6)))
  629. X      (or (= (% (+ (nth 2 a) ldiff) 2) 0) (setq ldiff (1+ ldiff)))
  630. X      (setq guess (math-make-float (math-isqrt-small
  631. X                    (math-scale-int (nth 1 a) (- ldiff)))
  632. X                       (/ (+ (nth 2 a) ldiff) 2)))))
  633. X    (math-sqrt-float-iter a guess))
  634. )
  635. X
  636. (defun math-sqrt-float-iter (a guess)   ; [F F F]
  637. X  (math-working "sqrt" guess)
  638. X  (let ((g2 (math-mul-float (math-add-float guess (math-div-float a guess))
  639. X                '(float 5 -1))))
  640. X     (if (math-nearly-equal-float g2 guess)
  641. X     g2
  642. X       (math-sqrt-float-iter a g2)))
  643. )
  644. X
  645. ;;; True if A and B differ only in the last digit of precision.  [P F F]
  646. (defun math-nearly-equal-float (a b)
  647. X  (let ((ediff (- (nth 2 a) (nth 2 b))))
  648. X    (cond ((= ediff 0)   ;; Expanded out for speed
  649. X       (setq ediff (math-add (Math-integer-neg (nth 1 a)) (nth 1 b)))
  650. X       (or (eq ediff 0)
  651. X           (and (not (consp ediff))
  652. X            (< ediff 10)
  653. X            (> ediff -10)
  654. X            (= (math-numdigs (nth 1 a)) calc-internal-prec))))
  655. X      ((= ediff 1)
  656. X       (setq ediff (math-add (Math-integer-neg (nth 1 b))
  657. X                 (math-scale-int (nth 1 a) 1)))
  658. X       (and (not (consp ediff))
  659. X        (< ediff 10)
  660. X        (> ediff -10)
  661. X        (= (math-numdigs (nth 1 b)) calc-internal-prec)))
  662. X      ((= ediff -1)
  663. X       (setq ediff (math-add (Math-integer-neg (nth 1 a))
  664. X                 (math-scale-int (nth 1 b) 1)))
  665. X       (and (not (consp ediff))
  666. X        (< ediff 10)
  667. X        (> ediff -10)
  668. X        (= (math-numdigs (nth 1 a)) calc-internal-prec)))))
  669. )
  670. X
  671. (defun math-nearly-equal (a b)   ;  [P N N] [Public]
  672. X  (setq a (math-float a))
  673. X  (setq b (math-float b))
  674. X  (if (eq (car a) 'polar) (setq a (math-complex a)))
  675. X  (if (eq (car b) 'polar) (setq b (math-complex b)))
  676. X  (if (eq (car a) 'cplx)
  677. X      (if (eq (car b) 'cplx)
  678. X      (and (or (math-nearly-equal-float (nth 1 a) (nth 1 b))
  679. X           (and (math-nearly-zerop-float (nth 1 a) (nth 2 a))
  680. X            (math-nearly-zerop-float (nth 1 b) (nth 2 b))))
  681. X           (or (math-nearly-equal-float (nth 2 a) (nth 2 b))
  682. X           (and (math-nearly-zerop-float (nth 2 a) (nth 1 a))
  683. X            (math-nearly-zerop-float (nth 2 b) (nth 1 b)))))
  684. X    (and (math-nearly-equal-float (nth 1 a) b)
  685. X         (math-nearly-zerop-float (nth 2 a) b)))
  686. X      (if (eq (car b) 'cplx)
  687. X      (and (math-nearly-equal-float a (nth 1 b))
  688. X           (math-nearly-zerop-float a (nth 2 b)))
  689. X    (math-nearly-equal-float a b)))
  690. )
  691. X
  692. ;;; True if A is nearly zero compared to B.  [P F F]
  693. (defun math-nearly-zerop-float (a b)
  694. X  (or (eq (nth 1 a) 0)
  695. X      (<= (+ (math-numdigs (nth 1 a)) (nth 2 a))
  696. X      (1+ (- (+ (math-numdigs (nth 1 b)) (nth 2 b)) calc-internal-prec))))
  697. )
  698. X
  699. (defun math-nearly-zerop (a b)   ; [P N R] [Public]
  700. X  (setq a (math-float a))
  701. X  (setq b (math-float b))
  702. X  (if (eq (car a) 'cplx)
  703. X      (and (math-nearly-zerop-float (nth 1 a) b)
  704. X       (math-nearly-zerop-float (nth 2 a) b))
  705. X    (if (eq (car a) 'polar)
  706. X    (math-nearly-zerop-float (nth 1 a) b)
  707. X      (math-nearly-zerop-float a b)))
  708. )
  709. X
  710. ;;; This implementation could be improved, accuracy-wise.
  711. (defun math-hypot (a b)
  712. X  (cond ((Math-zerop a) (math-abs b))
  713. X    ((Math-zerop b) (math-abs a))
  714. X    ((not (Math-scalarp a))
  715. X     (if (math-infinitep a)
  716. X         (if (math-infinitep b)
  717. X         (if (equal a b)
  718. X             a
  719. X           '(var nan var-nan))
  720. X           a)
  721. X       (calc-record-why 'scalarp a)
  722. X       (list 'calcFunc-hypot a b)))
  723. X    ((not (Math-scalarp b))
  724. X     (if (math-infinitep b)
  725. X         b
  726. X       (calc-record-why 'scalarp b)
  727. X       (list 'calcFunc-hypot a b)))
  728. X    ((and (Math-numberp a) (Math-numberp b))
  729. X     (math-with-extra-prec 1
  730. X       (math-sqrt (math-add (calcFunc-abssqr a) (calcFunc-abssqr b)))))
  731. X    ((eq (car-safe a) 'hms)
  732. X     (if (eq (car-safe b) 'hms)   ; this helps sdev's of hms forms
  733. X         (math-to-hms (math-hypot (math-from-hms a 'deg)
  734. X                      (math-from-hms b 'deg)))
  735. X       (math-to-hms (math-hypot (math-from-hms a 'deg) b))))
  736. X    ((eq (car-safe b) 'hms)
  737. X     (math-to-hms (math-hypot a (math-from-hms b 'deg))))
  738. X    (t nil))
  739. )
  740. (fset 'calcFunc-hypot (symbol-function 'math-hypot))
  741. X
  742. (defun calcFunc-sqr (x)
  743. X  (math-pow x 2)
  744. )
  745. X
  746. X
  747. X
  748. (defun math-nth-root (a n)
  749. X  (cond ((= n 2) (math-sqrt a))
  750. X    ((Math-zerop a) a)
  751. X    ((Math-negp a) nil)
  752. X    ((Math-integerp a)
  753. X     (let ((root (math-nth-root-integer a n)))
  754. X       (if (car root)
  755. X           (cdr root)
  756. X         (and (not calc-symbolic-mode)
  757. X          (math-nth-root-float (math-float a) n
  758. X                       (math-float (cdr root)))))))
  759. X    ((eq (car-safe a) 'frac)
  760. X     (let* ((num-root (math-nth-root-integer (nth 1 a) n))
  761. X        (den-root (math-nth-root-integer (nth 2 a) n)))
  762. X       (if (and (car num-root) (car den-root))
  763. X           (list 'frac (cdr num-root) (cdr den-root))
  764. X         (and (not calc-symbolic-mode)
  765. X          (math-nth-root-float
  766. X           (math-float a) n
  767. X           (math-div-float (math-float (cdr num-root))
  768. X                   (math-float (cdr den-root))))))))
  769. X    ((eq (car-safe a) 'float)
  770. X     (and (not calc-symbolic-mode)
  771. X          (math-nth-root-float a n)))
  772. X    ((eq (car-safe a) 'polar)
  773. X     (let ((root (math-nth-root (nth 1 a) n)))
  774. X       (and root (list 'polar root (math-div (nth 2 a) n)))))
  775. X    (t nil))
  776. )
  777. X
  778. (defun math-nth-root-float (a n &optional guess)
  779. X  (math-inexact-result)
  780. X  (math-with-extra-prec 1
  781. X    (let ((nf (math-float n))
  782. X      (nfm1 (math-float (1- n))))
  783. X      (math-nth-root-float-iter a (or guess
  784. X                      (math-make-float
  785. X                       1 (/ (+ (math-numdigs (nth 1 a))
  786. X                           (nth 2 a)
  787. X                           (/ n 2))
  788. X                        n))))))
  789. )
  790. X
  791. (defun math-nth-root-float-iter (a guess)   ; uses "n", "nf", "nfm1"
  792. X  (math-working "root" guess)
  793. X  (let ((g2 (math-div-float (math-add-float (math-mul nfm1 guess)
  794. X                        (math-div-float
  795. X                         a (math-ipow guess (1- n))))
  796. X                nf)))
  797. X    (if (math-nearly-equal-float g2 guess)
  798. X    g2
  799. X      (math-nth-root-float-iter a g2)))
  800. )
  801. X
  802. (defun math-nth-root-integer (a n &optional guess)   ; [I I S]
  803. X  (math-nth-root-int-iter a (or guess
  804. X                (math-scale-int 1 (/ (+ (math-numdigs a)
  805. X                            (1- n))
  806. X                             n))))
  807. )
  808. X
  809. (defun math-nth-root-int-iter (a guess)   ; uses "n"
  810. X  (math-working "root" guess)
  811. X  (let* ((q (math-idivmod a (math-ipow guess (1- n))))
  812. X     (s (math-add (car q) (math-mul (1- n) guess)))
  813. X     (g2 (math-idivmod s n)))
  814. X    (if (Math-natnum-lessp (car g2) guess)
  815. X    (math-nth-root-int-iter a (car g2))
  816. X      (cons (and (equal (car g2) guess)
  817. X         (eq (cdr q) 0)
  818. X         (eq (cdr g2) 0))
  819. X        guess)))
  820. )
  821. X
  822. (defun calcFunc-nroot (x n)
  823. X  (calcFunc-pow x (if (integerp n)
  824. X              (math-make-frac 1 n)
  825. X            (math-div 1 n)))
  826. )
  827. X
  828. X
  829. X
  830. X
  831. ;;;; Transcendental functions.
  832. X
  833. ;;; All of these functions are defined on the complex plane.
  834. ;;; (Branch cuts, etc. follow Steele's Common Lisp book.)
  835. X
  836. ;;; Most functions increase calc-internal-prec by 2 digits, then round
  837. ;;; down afterward.  "-raw" functions use the current precision, require
  838. ;;; their arguments to be in float (or complex float) format, and always
  839. ;;; work in radians (where applicable).
  840. X
  841. (defun math-to-radians (a)   ; [N N]
  842. X  (cond ((eq (car-safe a) 'hms)
  843. X     (math-from-hms a 'rad))
  844. X    ((memq calc-angle-mode '(deg hms))
  845. X     (math-mul a (math-pi-over-180)))
  846. X    (t a))
  847. )
  848. X
  849. (defun math-from-radians (a)   ; [N N]
  850. X  (cond ((eq calc-angle-mode 'deg)
  851. X     (if (math-constp a)
  852. X         (math-div a (math-pi-over-180))
  853. X       (list 'calcFunc-deg a)))
  854. X    ((eq calc-angle-mode 'hms)
  855. X     (math-to-hms a 'rad))
  856. X    (t a))
  857. )
  858. X
  859. (defun math-to-radians-2 (a)   ; [N N]
  860. X  (cond ((eq (car-safe a) 'hms)
  861. X     (math-from-hms a 'rad))
  862. X    ((memq calc-angle-mode '(deg hms))
  863. X     (if calc-symbolic-mode
  864. X         (math-div (math-mul a '(var pi var-pi)) 180)
  865. X       (math-mul a (math-pi-over-180))))
  866. X    (t a))
  867. )
  868. X
  869. (defun math-from-radians-2 (a)   ; [N N]
  870. X  (cond ((memq calc-angle-mode '(deg hms))
  871. X     (if calc-symbolic-mode
  872. X         (math-div (math-mul 180 a) '(var pi var-pi))
  873. X       (math-div a (math-pi-over-180))))
  874. X    (t a))
  875. )
  876. X
  877. X
  878. X
  879. ;;; Sine, cosine, and tangent.
  880. X
  881. (defun calcFunc-sin (x)   ; [N N] [Public]
  882. X  (cond ((and (integerp x)
  883. X          (if (eq calc-angle-mode 'deg)
  884. X          (= (% x 90) 0)
  885. X        (= x 0)))
  886. X     (aref [0 1 0 -1] (math-mod (/ x 90) 4)))
  887. X    ((Math-scalarp x)
  888. X     (math-with-extra-prec 2
  889. X       (math-sin-raw (math-to-radians (math-float x)))))
  890. X    ((eq (car x) 'sdev)
  891. X     (if (math-constp x)
  892. X         (math-with-extra-prec 2
  893. X           (let* ((xx (math-to-radians (math-float (nth 1 x))))
  894. X              (xs (math-to-radians (math-float (nth 2 x))))
  895. X              (sc (math-sin-cos-raw xx)))
  896. X         (math-make-sdev (car sc) (math-mul xs (cdr sc)))))
  897. X       (math-make-sdev (calcFunc-sin (nth 1 x))
  898. X               (math-mul (nth 2 x) (calcFunc-cos (nth 1 x))))))
  899. X    ((and (eq (car x) 'intv) (math-intv-constp x))
  900. X     (calcFunc-cos (math-sub x (math-quarter-circle nil))))
  901. X    ((equal x '(var nan var-nan))
  902. X     x)
  903. X    (t (calc-record-why 'scalarp x)
  904. X       (list 'calcFunc-sin x)))
  905. )
  906. X
  907. (defun calcFunc-cos (x)   ; [N N] [Public]
  908. X  (cond ((and (integerp x)
  909. X          (if (eq calc-angle-mode 'deg)
  910. X          (= (% x 90) 0)
  911. X        (= x 0)))
  912. X     (aref [1 0 -1 0] (math-mod (/ x 90) 4)))
  913. X    ((Math-scalarp x)
  914. X     (math-with-extra-prec 2
  915. X       (math-cos-raw (math-to-radians (math-float x)))))
  916. X    ((eq (car x) 'sdev)
  917. X     (if (math-constp x)
  918. X         (math-with-extra-prec 2
  919. X           (let* ((xx (math-to-radians (math-float (nth 1 x))))
  920. X              (xs (math-to-radians (math-float (nth 2 x))))
  921. X              (sc (math-sin-cos-raw xx)))
  922. X         (math-make-sdev (cdr sc) (math-mul xs (car sc)))))
  923. X       (math-make-sdev (calcFunc-cos (nth 1 x))
  924. X               (math-mul (nth 2 x) (calcFunc-sin (nth 1 x))))))
  925. X    ((and (eq (car x) 'intv) (math-intv-constp x))
  926. X     (math-with-extra-prec 2
  927. X       (let* ((xx (math-to-radians (math-float x)))
  928. X          (na (math-floor (math-div (nth 2 xx) (math-pi))))
  929. X          (nb (math-floor (math-div (nth 3 xx) (math-pi))))
  930. X          (span (math-sub nb na)))
  931. X         (if (memq span '(0 1))
  932. X         (let ((int (math-sort-intv (nth 1 x)
  933. X                        (math-cos-raw (nth 2 xx))
  934. X                        (math-cos-raw (nth 3 xx)))))
  935. X           (if (eq span 1)
  936. X               (if (math-evenp na)
  937. X               (math-make-intv (logior (nth 1 x) 2)
  938. X                       -1
  939. X                       (nth 3 int))
  940. X             (math-make-intv (logior (nth 1 x) 1)
  941. X                     (nth 2 int)
  942. X                     1))
  943. X             int))
  944. X           (list 'intv 3 -1 1)))))
  945. X    ((equal x '(var nan var-nan))
  946. X     x)
  947. X    (t (calc-record-why 'scalarp x)
  948. X       (list 'calcFunc-cos x)))
  949. )
  950. X
  951. (defun calcFunc-sincos (x)   ; [V N] [Public]
  952. X  (if (Math-scalarp x)
  953. X      (math-with-extra-prec 2
  954. X    (let ((sc (math-sin-cos-raw (math-to-radians (math-float x)))))
  955. X      (list 'vec (cdr sc) (car sc))))    ; the vector [cos, sin]
  956. X    (list 'vec (calcFunc-sin x) (calcFunc-cos x)))
  957. )
  958. X
  959. (defun calcFunc-tan (x)   ; [N N] [Public]
  960. X  (cond ((and (integerp x)
  961. X          (if (eq calc-angle-mode 'deg)
  962. X          (= (% x 180) 0)
  963. X        (= x 0)))
  964. X     0)
  965. X    ((Math-scalarp x)
  966. X     (math-with-extra-prec 2
  967. X       (math-tan-raw (math-to-radians (math-float x)))))
  968. X    ((eq (car x) 'sdev)
  969. X     (if (math-constp x)
  970. X         (math-with-extra-prec 2
  971. X           (let* ((xx (math-to-radians (math-float (nth 1 x))))
  972. X              (xs (math-to-radians (math-float (nth 2 x))))
  973. X              (sc (math-sin-cos-raw xx)))
  974. X         (if (and (math-zerop (cdr sc)) (not calc-infinite-mode))
  975. X             (progn
  976. X               (calc-record-why "*Division by zero")
  977. X               (list 'calcFunc-tan x))
  978. X           (math-make-sdev (math-div-float (car sc) (cdr sc))
  979. X                   (math-div-float xs (math-sqr (cdr sc)))))))
  980. X       (math-make-sdev (calcFunc-tan (nth 1 x))
  981. X               (math-div (nth 2 x)
  982. X                     (math-sqr (calcFunc-cos (nth 1 x)))))))
  983. X    ((and (eq (car x) 'intv) (math-intv-constp x))
  984. X     (or (math-with-extra-prec 2
  985. X           (let* ((xx (math-to-radians (math-float x)))
  986. X              (na (math-floor (math-div (math-sub (nth 2 xx)
  987. X                              (math-pi-over-2))
  988. X                        (math-pi))))
  989. X              (nb (math-floor (math-div (math-sub (nth 3 xx)
  990. X                              (math-pi-over-2))
  991. X                        (math-pi)))))
  992. X         (and (equal na nb)
  993. X              (math-sort-intv (nth 1 x)
  994. X                      (math-tan-raw (nth 2 xx))
  995. X                      (math-tan-raw (nth 3 xx))))))
  996. X         '(intv 3 (neg (var inf var-inf)) (var inf var-inf))))
  997. X    ((equal x '(var nan var-nan))
  998. X     x)
  999. X    (t (calc-record-why 'scalarp x)
  1000. X       (list 'calcFunc-tan x)))
  1001. )
  1002. X
  1003. (defun math-sin-raw (x)   ; [N N]
  1004. X  (cond ((eq (car x) 'cplx)
  1005. X     (let* ((expx (math-exp-raw (nth 2 x)))
  1006. X        (expmx (math-div-float '(float 1 0) expx))
  1007. X        (sc (math-sin-cos-raw (nth 1 x))))
  1008. X       (list 'cplx
  1009. X         (math-mul-float (car sc)
  1010. X                 (math-mul-float (math-add-float expx expmx)
  1011. X                         '(float 5 -1)))
  1012. X         (math-mul-float (cdr sc)
  1013. X                 (math-mul-float (math-sub-float expx expmx)
  1014. X                         '(float 5 -1))))))
  1015. X    ((eq (car x) 'polar)
  1016. X     (math-polar (math-sin-raw (math-complex x))))
  1017. X    ((Math-integer-negp (nth 1 x))
  1018. X     (math-neg-float (math-sin-raw (math-neg-float x))))
  1019. X    ((math-lessp-float '(float 7 0) x)  ; avoid inf loops due to roundoff
  1020. X     (math-sin-raw (math-mod x (math-two-pi))))
  1021. X    (t (math-sin-raw-2 x x)))
  1022. )
  1023. X
  1024. (defun math-cos-raw (x)   ; [N N]
  1025. X  (if (eq (car-safe x) 'polar)
  1026. X      (math-polar (math-cos-raw (math-complex x)))
  1027. X    (math-sin-raw (math-sub (math-pi-over-2) x)))
  1028. )
  1029. X
  1030. ;;; This could use a smarter method:  Reduce x as in math-sin-raw, then
  1031. ;;;   compute either sin(x) or cos(x), whichever is smaller, and compute
  1032. ;;;   the other using the identity sin(x)^2 + cos(x)^2 = 1.
  1033. (defun math-sin-cos-raw (x)   ; [F.F F]  (result is (sin x . cos x))
  1034. X  (cons (math-sin-raw x) (math-cos-raw x))
  1035. )
  1036. X
  1037. (defun math-tan-raw (x)   ; [N N]
  1038. X  (cond ((eq (car x) 'cplx)
  1039. X     (let* ((x (math-mul x '(float 2 0)))
  1040. X        (expx (math-exp-raw (nth 2 x)))
  1041. X        (expmx (math-div-float '(float 1 0) expx))
  1042. X        (sc (math-sin-cos-raw (nth 1 x)))
  1043. X        (d (math-add-float (cdr sc)
  1044. X                   (math-mul-float (math-add-float expx expmx)
  1045. X                           '(float 5 -1)))))
  1046. X       (and (not (eq (nth 1 d) 0))
  1047. X        (list 'cplx
  1048. X              (math-div-float (car sc) d)
  1049. X              (math-div-float (math-mul-float (math-sub-float expx
  1050. X                                      expmx)
  1051. X                              '(float 5 -1)) d)))))
  1052. X    ((eq (car x) 'polar)
  1053. X     (math-polar (math-tan-raw (math-complex x))))
  1054. X    (t
  1055. X     (let ((sc (math-sin-cos-raw x)))
  1056. X       (if (eq (nth 1 (cdr sc)) 0)
  1057. X           (math-div (car sc) 0)
  1058. X         (math-div-float (car sc) (cdr sc))))))
  1059. )
  1060. X
  1061. (defun math-sin-raw-2 (x orgx)   ; This avoids poss of inf recursion.  [F F]
  1062. X  (let ((xmpo2 (math-sub-float (math-pi-over-2) x)))
  1063. X    (cond ((Math-integer-negp (nth 1 xmpo2))
  1064. X       (math-neg-float (math-sin-raw-2 (math-sub-float x (math-pi))
  1065. X                       orgx)))
  1066. X      ((math-lessp-float (math-pi-over-4) x)
  1067. X       (math-cos-raw-2 xmpo2 orgx))
  1068. X      ((math-lessp-float x (math-neg (math-pi-over-4)))
  1069. X       (math-neg (math-cos-raw-2 (math-add (math-pi-over-2) x) orgx)))
  1070. X      ((math-nearly-zerop-float x orgx) '(float 0 0))
  1071. X      (calc-symbolic-mode (signal 'inexact-result nil))
  1072. X      (t (math-sin-series x 6 4 x (math-neg-float (math-sqr-float x))))))
  1073. )
  1074. X
  1075. (defun math-cos-raw-2 (x orgx)   ; [F F]
  1076. X  (cond ((math-nearly-zerop-float x orgx) '(float 1 0))
  1077. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1078. X    (t (let ((xnegsqr (math-neg-float (math-sqr-float x))))
  1079. X         (math-sin-series
  1080. X          (math-add-float '(float 1 0)
  1081. X                  (math-mul-float xnegsqr '(float 5 -1)))
  1082. X          24 5 xnegsqr xnegsqr))))
  1083. )
  1084. X
  1085. (defun math-sin-series (sum nfac n x xnegsqr)
  1086. X  (math-working "sin" sum)
  1087. X  (let* ((nextx (math-mul-float x xnegsqr))
  1088. X     (nextsum (math-add-float sum (math-div-float nextx
  1089. X                              (math-float nfac)))))
  1090. X    (if (math-nearly-equal-float sum nextsum)
  1091. X    sum
  1092. X      (math-sin-series nextsum (math-mul nfac (* n (1+ n)))
  1093. X               (+ n 2) nextx xnegsqr)))
  1094. )
  1095. X
  1096. X
  1097. ;;; Inverse sine, cosine, tangent.
  1098. X
  1099. (defun calcFunc-arcsin (x)   ; [N N] [Public]
  1100. X  (cond ((eq x 0) 0)
  1101. X    ((and (eq x 1) (eq calc-angle-mode 'deg)) 90)
  1102. X    ((and (eq x -1) (eq calc-angle-mode 'deg)) -90)
  1103. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1104. X    ((Math-numberp x)
  1105. X     (math-with-extra-prec 2
  1106. X       (math-from-radians (math-arcsin-raw (math-float x)))))
  1107. X    ((eq (car x) 'sdev)
  1108. X     (math-make-sdev (calcFunc-arcsin (nth 1 x))
  1109. X             (math-from-radians
  1110. X              (math-div (nth 2 x)
  1111. X                    (math-sqrt
  1112. X                     (math-sub 1 (math-sqr (nth 1 x))))))))
  1113. X    ((eq (car x) 'intv)
  1114. X     (math-sort-intv (nth 1 x)
  1115. X             (calcFunc-arcsin (nth 2 x))
  1116. X             (calcFunc-arcsin (nth 3 x))))
  1117. X    ((equal x '(var nan var-nan))
  1118. X     x)
  1119. X    (t (calc-record-why 'numberp x)
  1120. X       (list 'calcFunc-arcsin x)))
  1121. )
  1122. X
  1123. (defun calcFunc-arccos (x)   ; [N N] [Public]
  1124. X  (cond ((eq x 1) 0)
  1125. X    ((and (eq x 0) (eq calc-angle-mode 'deg)) 90)
  1126. X    ((and (eq x -1) (eq calc-angle-mode 'deg)) 180)
  1127. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1128. X    ((Math-numberp x)
  1129. X     (math-with-extra-prec 2
  1130. X       (math-from-radians (math-arccos-raw (math-float x)))))
  1131. X    ((eq (car x) 'sdev)
  1132. X     (math-make-sdev (calcFunc-arccos (nth 1 x))
  1133. X             (math-from-radians
  1134. X              (math-div (nth 2 x)
  1135. X                    (math-sqrt
  1136. X                     (math-sub 1 (math-sqr (nth 1 x))))))))
  1137. X    ((eq (car x) 'intv)
  1138. X     (math-sort-intv (nth 1 x)
  1139. X             (calcFunc-arccos (nth 2 x))
  1140. X             (calcFunc-arccos (nth 3 x))))
  1141. X    ((equal x '(var nan var-nan))
  1142. X     x)
  1143. X    (t (calc-record-why 'numberp x)
  1144. X       (list 'calcFunc-arccos x)))
  1145. )
  1146. X
  1147. (defun calcFunc-arctan (x)   ; [N N] [Public]
  1148. X  (cond ((eq x 0) 0)
  1149. X    ((and (eq x 1) (eq calc-angle-mode 'deg)) 45)
  1150. X    ((and (eq x -1) (eq calc-angle-mode 'deg)) -45)
  1151. X    ((Math-numberp x)
  1152. X     (math-with-extra-prec 2
  1153. X       (math-from-radians (math-arctan-raw (math-float x)))))
  1154. X    ((eq (car x) 'sdev)
  1155. X     (math-make-sdev (calcFunc-arctan (nth 1 x))
  1156. X             (math-from-radians
  1157. X              (math-div (nth 2 x)
  1158. X                    (math-add 1 (math-sqr (nth 1 x)))))))
  1159. X    ((eq (car x) 'intv)
  1160. X     (math-sort-intv (nth 1 x)
  1161. X             (calcFunc-arctan (nth 2 x))
  1162. X             (calcFunc-arctan (nth 3 x))))
  1163. X    ((equal x '(var inf var-inf))
  1164. X     (math-quarter-circle t))
  1165. X    ((equal x '(neg (var inf var-inf)))
  1166. X     (math-neg (math-quarter-circle t)))
  1167. X    ((equal x '(var nan var-nan))
  1168. X     x)
  1169. X    (t (calc-record-why 'numberp x)
  1170. X       (list 'calcFunc-arctan x)))
  1171. )
  1172. X
  1173. (defun math-arcsin-raw (x)   ; [N N]
  1174. X  (let ((a (math-sqrt-raw (math-sub '(float 1 0) (math-sqr x)))))
  1175. X    (if (or (memq (car x) '(cplx polar))
  1176. X        (memq (car a) '(cplx polar)))
  1177. X    (math-with-extra-prec 2   ; use extra precision for difficult case
  1178. X      (math-mul '(cplx 0 -1)
  1179. X            (math-ln-raw (math-add (math-mul '(cplx 0 1) x) a))))
  1180. X      (math-arctan2-raw x a)))
  1181. )
  1182. X
  1183. (defun math-arccos-raw (x)   ; [N N]
  1184. X  (math-sub (math-pi-over-2) (math-arcsin-raw x))
  1185. )
  1186. X
  1187. (defun math-arctan-raw (x)   ; [N N]
  1188. X  (cond ((memq (car x) '(cplx polar))
  1189. X     (math-with-extra-prec 2   ; extra-extra
  1190. X       (math-div (math-sub
  1191. X              (math-ln-raw (math-add 1 (math-mul '(cplx 0 1) x)))
  1192. X              (math-ln-raw (math-add 1 (math-mul '(cplx 0 -1) x))))
  1193. X             '(cplx 0 2))))
  1194. X    ((Math-integer-negp (nth 1 x))
  1195. X     (math-neg-float (math-arctan-raw (math-neg-float x))))
  1196. X    ((math-zerop x) x)
  1197. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1198. X    ((math-equal-int x 1) (math-pi-over-4))
  1199. X    ((math-equal-int x -1) (math-neg (math-pi-over-4)))
  1200. X    ((math-lessp-float '(float 414214 -6) x)  ; if x > sqrt(2) - 1, reduce
  1201. X     (if (math-lessp-float '(float 1 0) x)
  1202. X         (math-sub-float (math-mul-float (math-pi) '(float 5 -1))
  1203. X                 (math-arctan-raw (math-div-float '(float 1 0) x)))
  1204. X       (math-sub-float (math-mul-float (math-pi) '(float 25 -2))
  1205. X               (math-arctan-raw (math-div-float
  1206. X                         (math-sub-float '(float 1 0) x)
  1207. X                         (math-add-float '(float 1 0)
  1208. X                                 x))))))
  1209. X    (t (math-arctan-series x 3 x (math-neg-float (math-sqr-float x)))))
  1210. )
  1211. X
  1212. (defun math-arctan-series (sum n x xnegsqr)
  1213. X  (math-working "arctan" sum)
  1214. X  (let* ((nextx (math-mul-float x xnegsqr))
  1215. X     (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
  1216. X    (if (math-nearly-equal-float sum nextsum)
  1217. X    sum
  1218. X      (math-arctan-series nextsum (+ n 2) nextx xnegsqr)))
  1219. )
  1220. X
  1221. (defun calcFunc-arctan2 (y x)   ; [F R R] [Public]
  1222. X  (if (Math-anglep y)
  1223. X      (if (Math-anglep x)
  1224. X      (math-with-extra-prec 2
  1225. X        (math-from-radians (math-arctan2-raw (math-float y)
  1226. X                         (math-float x))))
  1227. X    (calc-record-why 'anglep x)
  1228. X    (list 'calcFunc-arctan2 y x))
  1229. X    (if (and (or (math-infinitep x) (math-anglep x))
  1230. X         (or (math-infinitep y) (math-anglep y)))
  1231. X    (progn
  1232. X      (if (math-posp x)
  1233. X          (setq x 1)
  1234. X        (if (math-negp x)
  1235. X        (setq x -1)
  1236. X          (or (math-zerop x)
  1237. X          (setq x nil))))
  1238. X      (if (math-posp y)
  1239. X          (setq y 1)
  1240. X        (if (math-negp y)
  1241. X        (setq y -1)
  1242. X          (or (math-zerop y)
  1243. X          (setq y nil))))
  1244. X      (if (and y x)
  1245. X          (calcFunc-arctan2 y x)
  1246. X        '(var nan var-nan)))
  1247. X      (calc-record-why 'anglep y)
  1248. X      (list 'calcFunc-arctan2 y x)))
  1249. )
  1250. X
  1251. (defun math-arctan2-raw (y x)   ; [F R R]
  1252. X  (cond ((math-zerop y)
  1253. X     (if (math-negp x) (math-pi)
  1254. X       (if (or (math-floatp x) (math-floatp y)) '(float 0 0) 0)))
  1255. X    ((math-zerop x)
  1256. X     (if (math-posp y)
  1257. X         (math-pi-over-2)
  1258. X       (math-neg (math-pi-over-2))))
  1259. X    ((math-posp x)
  1260. X     (math-arctan-raw (math-div-float y x)))
  1261. X    ((math-posp y)
  1262. X     (math-add-float (math-arctan-raw (math-div-float y x))
  1263. X             (math-pi)))
  1264. X    (t
  1265. X     (math-sub-float (math-arctan-raw (math-div-float y x))
  1266. X             (math-pi))))
  1267. )
  1268. X
  1269. (defun calcFunc-arcsincos (x)   ; [V N] [Public]
  1270. X  (if (and (Math-vectorp x)
  1271. X       (= (length x) 3))
  1272. X      (calcFunc-arctan2 (nth 2 x) (nth 1 x))
  1273. X    (math-reject-arg x "*Two-element vector expected"))
  1274. )
  1275. X
  1276. X
  1277. X
  1278. ;;; Exponential function.
  1279. X
  1280. (defun calcFunc-exp (x)   ; [N N] [Public]
  1281. X  (cond ((eq x 0) 1)
  1282. X    ((and (memq x '(1 -1)) calc-symbolic-mode)
  1283. X     (if (eq x 1) '(var e var-e) (math-div 1 '(var e var-e))))
  1284. X    ((Math-numberp x)
  1285. X     (math-with-extra-prec 2 (math-exp-raw (math-float x))))
  1286. X    ((eq (car-safe x) 'sdev)
  1287. X     (let ((ex (calcFunc-exp (nth 1 x))))
  1288. X       (math-make-sdev ex (math-mul (nth 2 x) ex))))
  1289. X    ((eq (car-safe x) 'intv)
  1290. X     (math-make-intv (nth 1 x) (calcFunc-exp (nth 2 x))
  1291. X             (calcFunc-exp (nth 3 x))))
  1292. X    ((equal x '(var inf var-inf))
  1293. X     x)
  1294. X    ((equal x '(neg (var inf var-inf)))
  1295. X     0)
  1296. X    ((equal x '(var nan var-nan))
  1297. X     x)
  1298. X    (t (calc-record-why 'numberp x)
  1299. X       (list 'calcFunc-exp x)))
  1300. )
  1301. X
  1302. (defun calcFunc-expm1 (x)   ; [N N] [Public]
  1303. X  (cond ((eq x 0) 0)
  1304. X    ((math-zerop x) '(float 0 0))
  1305. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1306. X    ((Math-numberp x)
  1307. X     (math-with-extra-prec 2
  1308. X       (let ((x (math-float x)))
  1309. X         (if (and (eq (car x) 'float)
  1310. X              (math-lessp-float x '(float 1 0))
  1311. X              (math-lessp-float '(float -1 0) x))
  1312. X         (math-exp-minus-1-raw x)
  1313. X           (math-add (math-exp-raw x) -1)))))
  1314. X    ((eq (car-safe x) 'sdev)
  1315. X     (if (math-constp x)
  1316. X         (let ((ex (calcFunc-expm1 (nth 1 x))))
  1317. X           (math-make-sdev ex (math-mul (nth 2 x) (math-add ex 1))))
  1318. X       (math-make-sdev (calcFunc-expm1 (nth 1 x))
  1319. X               (math-mul (nth 2 x) (calcFunc-exp (nth 1 x))))))
  1320. X    ((eq (car-safe x) 'intv)
  1321. X     (math-make-intv (nth 1 x)
  1322. X             (calcFunc-expm1 (nth 2 x))
  1323. X             (calcFunc-expm1 (nth 3 x))))
  1324. X    ((equal x '(var inf var-inf))
  1325. X     x)
  1326. X    ((equal x '(neg (var inf var-inf)))
  1327. X     -1)
  1328. X    ((equal x '(var nan var-nan))
  1329. X     x)
  1330. X    (t (calc-record-why 'numberp x)
  1331. X       (list 'calcFunc-expm1 x)))
  1332. )
  1333. X
  1334. (defun calcFunc-exp10 (x)   ; [N N] [Public]
  1335. X  (if (eq x 0)
  1336. X      1
  1337. X    (math-pow '(float 1 1) x))
  1338. )
  1339. X
  1340. (defun math-exp-raw (x)   ; [N N]
  1341. X  (cond ((math-zerop x) '(float 1 0))
  1342. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1343. X    ((eq (car x) 'cplx)
  1344. X     (let ((expx (math-exp-raw (nth 1 x)))
  1345. X           (sc (math-sin-cos-raw (nth 2 x))))
  1346. X       (list 'cplx
  1347. X         (math-mul-float expx (cdr sc))
  1348. X         (math-mul-float expx (car sc)))))
  1349. X    ((eq (car x) 'polar)
  1350. X     (let ((xc (math-complex x)))
  1351. X       (list 'polar
  1352. X         (math-exp-raw (nth 1 xc))
  1353. X         (math-from-radians (nth 2 xc)))))
  1354. X    ((or (math-lessp-float '(float 5 -1) x)
  1355. X         (math-lessp-float x '(float -5 -1)))
  1356. X     (if (math-lessp-float '(float 921035 1) x)
  1357. X         (math-overflow)
  1358. X       (if (math-lessp-float x '(float -921035 1))
  1359. X           (math-underflow)))
  1360. X     (let* ((two-x (math-mul-float x '(float 2 0)))
  1361. X        (hint (math-scale-int (nth 1 two-x) (nth 2 two-x)))
  1362. X        (hfrac (math-sub-float x (math-mul-float (math-float hint)
  1363. X                             '(float 5 -1)))))
  1364. X       (math-mul-float (math-ipow (math-sqrt-e) hint)
  1365. X               (math-add-float '(float 1 0)
  1366. X                       (math-exp-minus-1-raw hfrac)))))
  1367. X    (t (math-add-float '(float 1 0) (math-exp-minus-1-raw x))))
  1368. )
  1369. X
  1370. (defun math-exp-minus-1-raw (x)   ; [F F]
  1371. X  (math-exp-series x 2 3 x x)
  1372. )
  1373. X
  1374. (defun math-exp-series (sum nfac n xpow x)
  1375. X  (math-working "exp" sum)
  1376. X  (let* ((nextx (math-mul-float xpow x))
  1377. X     (nextsum (math-add-float sum (math-div-float nextx
  1378. X                              (math-float nfac)))))
  1379. X    (if (math-nearly-equal-float sum nextsum)
  1380. X    sum
  1381. X      (math-exp-series nextsum (math-mul nfac n) (1+ n) nextx x)))
  1382. )
  1383. X
  1384. X
  1385. X
  1386. ;;; Logarithms.
  1387. X
  1388. (defun calcFunc-ln (x)   ; [N N] [Public]
  1389. X  (cond ((math-zerop x)
  1390. X     (if calc-infinite-mode
  1391. X         '(neg (var inf var-inf))
  1392. X       (math-reject-arg x "*Logarithm of zero")))
  1393. X    ((eq x 1) 0)
  1394. X    ((Math-numberp x)
  1395. X     (math-with-extra-prec 2 (math-ln-raw (math-float x))))
  1396. X    ((eq (car-safe x) 'sdev)
  1397. X     (math-make-sdev (calcFunc-ln (nth 1 x))
  1398. X             (math-div (nth 2 x) (nth 1 x))))
  1399. X    ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
  1400. X                      (Math-zerop (nth 2 x))
  1401. X                      (not (math-intv-constp x))))
  1402. X     (let ((calc-infinite-mode t))
  1403. X       (math-make-intv (nth 1 x) (calcFunc-ln (nth 2 x))
  1404. X               (calcFunc-ln (nth 3 x)))))
  1405. X    ((equal x '(var e var-e))
  1406. X     1)
  1407. X    ((and (eq (car-safe x) '^)
  1408. X          (equal (nth 1 x) '(var e var-e))
  1409. X          (math-known-realp (nth 2 x)))
  1410. X     (nth 2 x))
  1411. X    ((math-infinitep x)
  1412. X     (if (equal x '(var nan var-nan))
  1413. X         x
  1414. X       '(var inf var-inf)))
  1415. X    (t (calc-record-why 'numberp x)
  1416. X       (list 'calcFunc-ln x)))
  1417. )
  1418. X
  1419. (defun calcFunc-log10 (x)   ; [N N] [Public]
  1420. X  (cond ((math-equal-int x 1)
  1421. X     (if (math-floatp x) '(float 0 0) 0))
  1422. X    ((and (Math-integerp x)
  1423. X          (math-posp x)
  1424. X          (let ((res (math-integer-log x 10)))
  1425. X        (and (car res)
  1426. X             (setq x (cdr res)))))
  1427. X     x)
  1428. X    ((and (eq (car-safe x) 'frac)
  1429. X          (eq (nth 1 x) 1)
  1430. X          (let ((res (math-integer-log (nth 2 x) 10)))
  1431. X        (and (car res)
  1432. X             (setq x (- (cdr res))))))
  1433. X     x)
  1434. X    ((math-zerop x)
  1435. X     (if calc-infinite-mode
  1436. X         '(neg (var inf var-inf))
  1437. X       (math-reject-arg x "*Logarithm of zero")))
  1438. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1439. X    ((Math-numberp x)
  1440. X     (math-with-extra-prec 2
  1441. X       (let ((xf (math-float x)))
  1442. X         (if (eq (nth 1 xf) 0)
  1443. X         (math-reject-arg x "*Logarithm of zero"))
  1444. X         (if (Math-integer-posp (nth 1 xf))
  1445. X         (if (eq (nth 1 xf) 1)    ; log10(1*10^n) = n
  1446. X             (math-float (nth 2 xf))
  1447. X           (let ((xdigs (1- (math-numdigs (nth 1 xf)))))
  1448. X             (math-add-float
  1449. X              (math-div-float (math-ln-raw-2
  1450. X                       (list 'float (nth 1 xf) (- xdigs)))
  1451. X                      (math-ln-10))
  1452. X              (math-float (+ (nth 2 xf) xdigs)))))
  1453. X           (math-div (calcFunc-ln xf) (math-ln-10))))))
  1454. X    ((eq (car-safe x) 'sdev)
  1455. X     (math-make-sdev (calcFunc-log10 (nth 1 x))
  1456. X             (math-div (nth 2 x)
  1457. X                   (math-mul (nth 1 x) (math-ln-10)))))
  1458. X    ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
  1459. X                      (not (math-intv-constp x))))
  1460. X     (math-make-intv (nth 1 x)
  1461. X             (calcFunc-log10 (nth 2 x))
  1462. X             (calcFunc-log10 (nth 3 x))))
  1463. X    ((math-infinitep x)
  1464. X     (if (equal x '(var nan var-nan))
  1465. X         x
  1466. X       '(var inf var-inf)))
  1467. X    (t (calc-record-why 'numberp x)
  1468. X       (list 'calcFunc-log10 x)))
  1469. )
  1470. X
  1471. (defun calcFunc-log (x &optional b)   ; [N N N] [Public]
  1472. X  (cond ((or (null b) (equal b '(var e var-e)))
  1473. X     (math-normalize (list 'calcFunc-ln x)))
  1474. X    ((or (eq b 10) (equal b '(float 1 1)))
  1475. X     (math-normalize (list 'calcFunc-log10 x)))
  1476. X    ((math-zerop x)
  1477. X     (if calc-infinite-mode
  1478. X         (math-div (calcFunc-ln x) (calcFunc-ln b))
  1479. X       (math-reject-arg x "*Logarithm of zero")))
  1480. X    ((math-zerop b)
  1481. X     (if calc-infinite-mode
  1482. X         (math-div (calcFunc-ln x) (calcFunc-ln b))
  1483. X       (math-reject-arg b "*Logarithm of zero")))
  1484. X    ((math-equal-int b 1)
  1485. X     (if calc-infinite-mode
  1486. X         (math-div (calcFunc-ln x) 0)
  1487. X       (math-reject-arg b "*Logarithm base one")))
  1488. X    ((math-equal-int x 1)
  1489. X     (if (or (math-floatp a) (math-floatp b)) '(float 0 0) 0))
  1490. X    ((and (Math-ratp x) (Math-ratp b)
  1491. X          (math-posp x) (math-posp b)
  1492. X          (let* ((sign 1) (inv nil)
  1493. X             (xx (if (Math-lessp 1 x)
  1494. X                 x
  1495. X               (setq sign -1)
  1496. X               (math-div 1 x)))
  1497. X             (bb (if (Math-lessp 1 b)
  1498. X                 b
  1499. X               (setq sign (- sign))
  1500. X               (math-div 1 b)))
  1501. X             (res (if (Math-lessp xx bb)
  1502. X                  (setq inv (math-integer-log bb xx))
  1503. X                (math-integer-log xx bb))))
  1504. X        (and (car res)
  1505. X             (setq x (if inv
  1506. X                 (math-div 1 (* sign (cdr res)))
  1507. X                   (* sign (cdr res)))))))
  1508. X     x)
  1509. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1510. X    ((and (Math-numberp x) (Math-numberp b))
  1511. X     (math-with-extra-prec 2
  1512. X       (math-div (math-ln-raw (math-float x))
  1513. X             (math-log-base-raw b))))
  1514. X    ((and (eq (car-safe x) 'sdev)
  1515. X          (Math-numberp b))
  1516. X     (math-make-sdev (calcFunc-log (nth 1 x) b)
  1517. X             (math-div (nth 2 x)
  1518. X                   (math-mul (nth 1 x)
  1519. X                         (math-log-base-raw b)))))
  1520. X    ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
  1521. X                      (not (math-intv-constp x)))
  1522. X          (math-realp b))
  1523. X     (math-make-intv (nth 1 x)
  1524. X             (calcFunc-log (nth 2 x) b)
  1525. X             (calcFunc-log (nth 3 x) b)))
  1526. X    ((or (eq (car-safe x) 'intv) (eq (car-safe b) 'intv))
  1527. X     (math-div (calcFunc-ln x) (calcFunc-ln b)))
  1528. X    ((or (math-infinitep x)
  1529. X         (math-infinitep b))
  1530. X     (math-div (calcFunc-ln x) (calcFunc-ln b)))
  1531. X    (t (if (Math-numberp b)
  1532. X           (calc-record-why 'numberp x)
  1533. X         (calc-record-why 'numberp b))
  1534. X       (list 'calcFunc-log x b)))
  1535. )
  1536. X
  1537. (defun calcFunc-alog (x &optional b)
  1538. X  (cond ((or (null b) (equal b '(var e var-e)))
  1539. X     (math-normalize (list 'calcFunc-exp x)))
  1540. X    (t (math-pow b x)))
  1541. )
  1542. X
  1543. (defun calcFunc-ilog (x b)
  1544. X  (if (and (math-natnump x) (not (eq x 0))
  1545. X       (math-natnump b) (not (eq b 0)))
  1546. X      (if (eq b 1)
  1547. X      (math-reject-arg x "*Logarithm base one")
  1548. X    (if (Math-natnum-lessp x b)
  1549. X        0
  1550. X      (cdr (math-integer-log x b))))
  1551. X    (math-floor (calcFunc-log x b)))
  1552. )
  1553. X
  1554. (defun math-integer-log (x b)
  1555. X  (let ((pows (list b))
  1556. X    (pow (math-sqr b))
  1557. X    next
  1558. X    sum n)
  1559. X    (while (not (Math-lessp x pow))
  1560. X      (setq pows (cons pow pows)
  1561. X        pow (math-sqr pow)))
  1562. X    (setq n (lsh 1 (1- (length pows)))
  1563. X      sum n
  1564. X      pow (car pows))
  1565. X    (while (and (setq pows (cdr pows))
  1566. X        (Math-lessp pow x))
  1567. X      (setq n (/ n 2)
  1568. X        next (math-mul pow (car pows)))
  1569. X      (or (Math-lessp x next)
  1570. X      (setq pow next
  1571. X        sum (+ sum n))))
  1572. X    (cons (equal pow x) sum))
  1573. )
  1574. X
  1575. X
  1576. (defun math-log-base-raw (b)   ; [N N]
  1577. X  (if (not (and (equal (car math-log-base-cache) b)
  1578. X        (eq (nth 1 math-log-base-cache) calc-internal-prec)))
  1579. X      (setq math-log-base-cache (list b calc-internal-prec
  1580. X                      (math-ln-raw (math-float b)))))
  1581. X  (nth 2 math-log-base-cache)
  1582. )
  1583. (setq math-log-base-cache nil)
  1584. X
  1585. (defun calcFunc-lnp1 (x)   ; [N N] [Public]
  1586. X  (cond ((Math-equal-int x -1)
  1587. X     (if calc-infinite-mode
  1588. X         '(neg (var inf var-inf))
  1589. X       (math-reject-arg x "*Logarithm of zero")))
  1590. X    ((eq x 0) 0)
  1591. X    ((math-zerop x) '(float 0 0))
  1592. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1593. X    ((Math-numberp x)
  1594. X     (math-with-extra-prec 2
  1595. X       (let ((x (math-float x)))
  1596. X         (if (and (eq (car x) 'float)
  1597. X              (math-lessp-float x '(float 5 -1))
  1598. X              (math-lessp-float '(float -5 -1) x))
  1599. X         (math-ln-plus-1-raw x)
  1600. X           (math-ln-raw (math-add-float x '(float 1 0)))))))
  1601. X    ((eq (car-safe x) 'sdev)
  1602. X     (math-make-sdev (calcFunc-lnp1 (nth 1 x))
  1603. X             (math-div (nth 2 x) (math-add (nth 1 x) 1))))
  1604. X    ((and (eq (car-safe x) 'intv) (or (Math-posp (nth 2 x))
  1605. X                      (not (math-intv-constp x))))
  1606. X     (math-make-intv (nth 1 x)
  1607. X             (calcFunc-lnp1 (nth 2 x))
  1608. X             (calcFunc-lnp1 (nth 3 x))))
  1609. X    ((math-infinitep x)
  1610. X     (if (equal x '(var nan var-nan))
  1611. X         x
  1612. X       '(var inf var-inf)))
  1613. X    (t (calc-record-why 'numberp x)
  1614. X       (list 'calcFunc-lnp1 x)))
  1615. )
  1616. X
  1617. (defun math-ln-raw (x)    ; [N N] --- must be float format!
  1618. X  (cond ((eq (car-safe x) 'cplx)
  1619. X     (list 'cplx
  1620. X           (math-mul-float (math-ln-raw
  1621. X                (math-add-float (math-sqr-float (nth 1 x))
  1622. X                        (math-sqr-float (nth 2 x))))
  1623. X                   '(float 5 -1))
  1624. X           (math-arctan2-raw (nth 2 x) (nth 1 x))))
  1625. X    ((eq (car x) 'polar)
  1626. X     (math-polar (list 'cplx
  1627. X               (math-ln-raw (nth 1 x))
  1628. X               (math-to-radians (nth 2 x)))))
  1629. X    ((Math-equal-int x 1)
  1630. X     '(float 0 0))
  1631. X    (calc-symbolic-mode (signal 'inexact-result nil))
  1632. X    ((math-posp (nth 1 x))    ; positive and real
  1633. X     (let ((xdigs (1- (math-numdigs (nth 1 x)))))
  1634. X       (math-add-float (math-ln-raw-2 (list 'float (nth 1 x) (- xdigs)))
  1635. X               (math-mul-float (math-float (+ (nth 2 x) xdigs))
  1636. X                       (math-ln-10)))))
  1637. X    ((math-zerop x)
  1638. X     (math-reject-arg x "*Logarithm of zero"))
  1639. X    ((eq calc-complex-mode 'polar)    ; negative and real
  1640. X     (math-polar
  1641. X      (list 'cplx   ; negative and real
  1642. X        (math-ln-raw (math-neg-float x))
  1643. X        (math-pi))))
  1644. X    (t (list 'cplx   ; negative and real
  1645. X         (math-ln-raw (math-neg-float x))
  1646. X         (math-pi))))
  1647. )
  1648. X
  1649. (defun math-ln-raw-2 (x)    ; [F F]
  1650. X  (cond ((math-lessp-float '(float 14 -1) x)
  1651. X     (math-add-float (math-ln-raw-2 (math-mul-float x '(float 5 -1)))
  1652. X             (math-ln-2)))
  1653. X    (t    ; now .7 < x <= 1.4
  1654. X     (math-ln-raw-3 (math-div-float (math-sub-float x '(float 1 0))
  1655. X                    (math-add-float x '(float 1 0))))))
  1656. )
  1657. X
  1658. (defun math-ln-raw-3 (x)   ; [F F]
  1659. X  (math-mul-float (math-ln-raw-series x 3 x (math-sqr-float x))
  1660. X          '(float 2 0))
  1661. )
  1662. X
  1663. ;;; Compute ln((1+x)/(1-x))
  1664. (defun math-ln-raw-series (sum n x xsqr)
  1665. X  (math-working "log" sum)
  1666. X  (let* ((nextx (math-mul-float x xsqr))
  1667. X     (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
  1668. X    (if (math-nearly-equal-float sum nextsum)
  1669. X    sum
  1670. X      (math-ln-raw-series nextsum (+ n 2) nextx xsqr)))
  1671. )
  1672. X
  1673. (defun math-ln-plus-1-raw (x)
  1674. X  (math-lnp1-series x 2 x (math-neg x))
  1675. )
  1676. X
  1677. (defun math-lnp1-series (sum n xpow x)
  1678. X  (math-working "lnp1" sum)
  1679. X  (let* ((nextx (math-mul-float xpow x))
  1680. X     (nextsum (math-add-float sum (math-div-float nextx (math-float n)))))
  1681. X    (if (math-nearly-equal-float sum nextsum)
  1682. X    sum
  1683. X      (math-lnp1-series nextsum (1+ n) nextx x)))
  1684. )
  1685. X
  1686. (math-defcache math-ln-10 (float (bigpos 018 684 045 994 092 585 302 2) -21)
  1687. X  (math-ln-raw-2 '(float 1 1)))
  1688. X
  1689. (math-defcache math-ln-2 (float (bigpos 417 309 945 559 180 147 693) -21)
  1690. X  (math-ln-raw-3 (math-float '(frac 1 3))))
  1691. X
  1692. X
  1693. X
  1694. ;;; Hyperbolic functions.
  1695. X
  1696. (defun calcFunc-sinh (x)   ; [N N] [Public]
  1697. X  (cond ((eq x 0) 0)
  1698. X    (math-expand-formulas
  1699. X     (math-normalize
  1700. X      (list '/ (list '- (list 'calcFunc-exp x)
  1701. X             (list 'calcFunc-exp (list 'neg x))) 2)))
  1702. X    ((Math-numberp x)
  1703. X     (if calc-symbolic-mode (signal 'inexact-result nil))
  1704. X     (math-with-extra-prec 2
  1705. X       (let ((expx (math-exp-raw (math-float x))))
  1706. X         (math-mul (math-add expx (math-div -1 expx)) '(float 5 -1)))))
  1707. X    ((eq (car-safe x) 'sdev)
  1708. X     (math-make-sdev (calcFunc-sinh (nth 1 x))
  1709. X             (math-mul (nth 2 x) (calcFunc-cosh (nth 1 x)))))
  1710. X    ((eq (car x) 'intv)
  1711. X     (math-sort-intv (nth 1 x)
  1712. X             (calcFunc-sinh (nth 2 x))
  1713. X             (calcFunc-sinh (nth 3 x))))
  1714. X    ((or (equal x '(var inf var-inf))
  1715. X         (equal x '(neg (var inf var-inf)))
  1716. X         (equal x '(var nan var-nan)))
  1717. X     x)
  1718. X    (t (calc-record-why 'numberp x)
  1719. X       (list 'calcFunc-sinh x)))
  1720. )
  1721. (put 'calcFunc-sinh 'math-expandable t)
  1722. X
  1723. (defun calcFunc-cosh (x)   ; [N N] [Public]
  1724. X  (cond ((eq x 0) 1)
  1725. X    (math-expand-formulas
  1726. X     (math-normalize
  1727. X      (list '/ (list '+ (list 'calcFunc-exp x)
  1728. X             (list 'calcFunc-exp (list 'neg x))) 2)))
  1729. X    ((Math-numberp x)
  1730. X     (if calc-symbolic-mode (signal 'inexact-result nil))
  1731. X     (math-with-extra-prec 2
  1732. X       (let ((expx (math-exp-raw (math-float x))))
  1733. X         (math-mul (math-add expx (math-div 1 expx)) '(float 5 -1)))))
  1734. X    ((eq (car-safe x) 'sdev)
  1735. X     (math-make-sdev (calcFunc-cosh (nth 1 x))
  1736. X             (math-mul (nth 2 x)
  1737. X                   (calcFunc-sinh (nth 1 x)))))
  1738. X    ((and (eq (car x) 'intv) (math-intv-constp x))
  1739. X     (setq x (math-abs x))
  1740. X     (math-sort-intv (nth 1 x)
  1741. X             (calcFunc-cosh (nth 2 x))
  1742. X             (calcFunc-cosh (nth 3 x))))
  1743. X    ((or (equal x '(var inf var-inf))
  1744. X         (equal x '(neg (var inf var-inf)))
  1745. X         (equal x '(var nan var-nan)))
  1746. X     (math-abs x))
  1747. X    (t (calc-record-why 'numberp x)
  1748. X       (list 'calcFunc-cosh x)))
  1749. )
  1750. (put 'calcFunc-cosh 'math-expandable t)
  1751. X
  1752. (defun calcFunc-tanh (x)   ; [N N] [Public]
  1753. X  (cond ((eq x 0) 0)
  1754. X    (math-expand-formulas
  1755. X     (math-normalize
  1756. X      (let ((expx (list 'calcFunc-exp x))
  1757. X        (expmx (list 'calcFunc-exp (list 'neg x))))
  1758. X        (math-normalize
  1759. X         (list '/ (list '- expx expmx) (list '+ expx expmx))))))
  1760. X    ((Math-numberp x)
  1761. X     (if calc-symbolic-mode (signal 'inexact-result nil))
  1762. X     (math-with-extra-prec 2
  1763. X       (let* ((expx (calcFunc-exp (math-float x)))
  1764. X          (expmx (math-div 1 expx)))
  1765. X         (math-div (math-sub expx expmx)
  1766. X               (math-add expx expmx)))))
  1767. X    ((eq (car-safe x) 'sdev)
  1768. X     (math-make-sdev (calcFunc-tanh (nth 1 x))
  1769. X             (math-div (nth 2 x)
  1770. X                   (math-sqr (calcFunc-cosh (nth 1 x))))))
  1771. X    ((eq (car x) 'intv)
  1772. X     (math-sort-intv (nth 1 x)
  1773. X             (calcFunc-tanh (nth 2 x))
  1774. X             (calcFunc-tanh (nth 3 x))))
  1775. X    ((equal x '(var inf var-inf))
  1776. X     1)
  1777. X    ((equal x '(neg (var inf var-inf)))
  1778. X     -1)
  1779. X    ((equal x '(var nan var-nan))
  1780. X     x)
  1781. X    (t (calc-record-why 'numberp x)
  1782. X       (list 'calcFunc-tanh x)))
  1783. )
  1784. (put 'calcFunc-tanh 'math-expandable t)
  1785. X
  1786. (defun calcFunc-arcsinh (x)   ; [N N] [Public]
  1787. X  (cond ((eq x 0) 0)
  1788. X    (math-expand-formulas
  1789. X     (math-normalize
  1790. X      (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
  1791. X                          (list '+ (list '^ x 2) 1))))))
  1792. X    ((Math-numberp x)
  1793. X     (if calc-symbolic-mode (signal 'inexact-result nil))
  1794. X     (math-with-extra-prec 2
  1795. X       (math-ln-raw (math-add x (math-sqrt-raw (math-add (math-sqr x)
  1796. X                                 '(float 1 0)))))))
  1797. X    ((eq (car-safe x) 'sdev)
  1798. X     (math-make-sdev (calcFunc-arcsinh (nth 1 x))
  1799. X             (math-div (nth 2 x)
  1800. X                   (math-sqrt
  1801. X                    (math-add (math-sqr (nth 1 x)) 1)))))
  1802. X    ((eq (car x) 'intv)
  1803. X     (math-sort-intv (nth 1 x)
  1804. X             (calcFunc-arcsinh (nth 2 x))
  1805. X             (calcFunc-arcsinh (nth 3 x))))
  1806. X    ((or (equal x '(var inf var-inf))
  1807. X         (equal x '(neg (var inf var-inf)))
  1808. X         (equal x '(var nan var-nan)))
  1809. X     x)
  1810. X    (t (calc-record-why 'numberp x)
  1811. X       (list 'calcFunc-arcsinh x)))
  1812. )
  1813. (put 'calcFunc-arcsinh 'math-expandable t)
  1814. X
  1815. (defun calcFunc-arccosh (x)   ; [N N] [Public]
  1816. X  (cond ((eq x 1) 0)
  1817. X    ((and (eq x -1) calc-symbolic-mode)
  1818. X     '(var pi var-pi))
  1819. X    ((and (eq x 0) calc-symbolic-mode)
  1820. X     (math-div (math-mul '(var pi var-pi) '(var i var-i)) 2))
  1821. X    (math-expand-formulas
  1822. X     (math-normalize
  1823. X      (list 'calcFunc-ln (list '+ x (list 'calcFunc-sqrt
  1824. X                          (list '- (list '^ x 2) 1))))))
  1825. X    ((Math-numberp x)
  1826. X     (if calc-symbolic-mode (signal 'inexact-result nil))
  1827. X     (if (Math-equal-int x -1)
  1828. X         (math-imaginary (math-pi))
  1829. X       (math-with-extra-prec 2
  1830. X         (if (or t    ; need to do this even in the real case!
  1831. X             (memq (car-safe x) '(cplx polar)))
  1832. X         (let ((xp1 (math-add 1 x)))  ; this gets the branch cuts right
  1833. X           (math-ln-raw
  1834. X            (math-add x (math-mul xp1
  1835. X                      (math-sqrt-raw
  1836. X                       (math-div (math-sub
  1837. X                              x
  1838. X                              '(float 1 0))
  1839. X                             xp1))))))
  1840. X           (math-ln-raw
  1841. X        (math-add x (math-sqrt-raw (math-add (math-sqr x)
  1842. X                             '(float -1 0)))))))))
  1843. X    ((eq (car-safe x) 'sdev)
  1844. X     (math-make-sdev (calcFunc-arccosh (nth 1 x))
  1845. X             (math-div (nth 2 x)
  1846. X                   (math-sqrt
  1847. X                    (math-add (math-sqr (nth 1 x)) -1)))))
  1848. X    ((eq (car x) 'intv)
  1849. X     (math-sort-intv (nth 1 x)
  1850. X             (calcFunc-arccosh (nth 2 x))
  1851. X             (calcFunc-arccosh (nth 3 x))))
  1852. X    ((or (equal x '(var inf var-inf))
  1853. X         (equal x '(neg (var inf var-inf)))
  1854. X         (equal x '(var nan var-nan)))
  1855. SHAR_EOF
  1856. true || echo 'restore of calc-math.el failed'
  1857. fi
  1858. echo 'End of  part 21'
  1859. echo 'File calc-math.el is continued in part 22'
  1860. echo 22 > _shar_seq_.tmp
  1861. exit 0
  1862. exit 0 # Just in case...
  1863. -- 
  1864. Kent Landfield                   INTERNET: kent@sparky.IMD.Sterling.COM
  1865. Sterling Software, IMD           UUCP:     uunet!sparky!kent
  1866. Phone:    (402) 291-8300         FAX:      (402) 291-4362
  1867. Please send comp.sources.misc-related mail to kent@uunet.uu.net.
  1868.